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Abstract. We use a 3D general circulation model to compare the primitive Martian 
hydrological cycle in “warm and wet” and “cold and icy” scenarios. In the “warm and 
wet” scenario, an anomalously high solar flux or intense greenhouse warming artificially 
added to the climate model are required to maintain warm conditions and an ice-free 
northern ocean. Precipitation shows strong surface variations, with high rates around Hel¬ 
las basin and west of Tlrarsis but low rates around Margaritifer Sinus (where the ob¬ 
served valley network drainage density is nonetheless high). In the “cold and icy” sce¬ 
nario, snow migration is a function of both obliquity and surface pressure, and limited 
episodic melting is possible through combinations of seasonal, volcanic and impact forc¬ 
ing. At surface pressures above those required to avoid atmospheric collapse (~ 0.5 bar) 
and moderate to high obliquity, snow is transported to the equatorial highland regions 
where the concentration of valley networks is highest. Snow accumulation in the Aeo- 
lis quadrangle is high, indicating an ice-free northern ocean is not required to supply wa¬ 
ter to Gale crater. At lower surface pressures and obliquities, both H 2 0 and C0 2 are 
trapped as ice at the poles and the equatorial regions become extremely dry. The val¬ 
ley network distribution is positively correlated with snow accumulation produced by the 
“cold and icy” simulation at 41.8° obliquity but uncorrelated with precipitation produced 
by the “warm and wet” simulation. Because our simulations make specific predictions 
for precipitation patterns under different climate scenarios, they motivate future targeted 
geological studies. 


1. Background 

Despite decades of research, deciphering the nature of 
Mars’ early climate remains a huge challenge. Although 
Mars receives only 43% of the solar flux incident on Earth, 
and the Sun’s luminosity was likely 20-30% lower 3-4 Ga, 
there is extensive evidence for aqueous alteration on Mars’ 
late Noachian and early Hesperian terrain. This evidence in¬ 
cludes dendritic valley networks (VNs) that are distributed 
widely across low to mid latitudes [Carr, 1996; Mangold 
et al., 2004; Hynek et al., 2010], open-basin lakes [Fas- 
sett and Head, 2008b], in-situ observations of conglomer¬ 
ates [Williams et al., 2013], and spectroscopic observations 
of phyllosilicate and sulphate minerals [Bibring et al., 2006; 
Mustard et al., 2008; Ehlmann et al., 2011]. 
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All these features indicate a pervasive influence of liq¬ 
uid water on the early Martian surface. Nonetheless, key 
uncertainties in the nature of the early Martian surface en¬ 
vironment remain. These include the intensity and dura¬ 
tion of warming episodes, and the extent to which the total 
surface H a O and C0 2 inventories were greater than today. 
Broadly speaking, proposed solutions to the problem can 
be divided into those that invoke long-term warm, wet con¬ 
ditions (e.g., Pollack et al. [1987]; Craddock and Howard 
[2002]), and those that assume the planet was mainly frozen, 
with aquifer discharge or episodic / seasonal melting of snow 
and ice deposits providing the necessary liquid water for flu¬ 
vial erosion (e.g., Squyres and Kasting [1994]; Toon et al. 
[2010]; Wordsworth et al. [2013]). 

Previously, we have shown that long-term warm, wet con¬ 
ditions on early Mars are not achieved even when the effects 
of clouds, dust and C0 2 collision-induced absorption (CIA) 
are taken into account in 3D models [Wordsworth et al., 
2013; Forget et al., 2013], because of a combination of the 
faint young Sun effect [Sagan and Mullen, 1972] and the con¬ 
densation of C0 2 at high pressures [Kasting, 1991]. Instead, 
we have proposed that the migration of ice towards high 
altitude regions due to adiabatic surface cooling (the ‘icy 
highlands’ effect) may have allowed the continual recharge 
of valley network sources after transient melting events in a 
predominately cold climate [Wordsworth et al., 2013]. 

2. Motivation 

The spatial distribution of the ancient fluvial features on 
Mars is strongly heterogeneous. In some cases, such as the 
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northern lowlands, the absence of surface evidence for liq¬ 
uid water is clearly a result of coverage by volcanic mate¬ 
rial and/or impact debris in the Hesperian and Amazonian 
epochs [Tanaka, 1986; Head et al., 2002], However, even 
in the Noachian highlands, the distribution of VNs [Hynek 
et al., 2010], open-basin lakes [Fassett and Head, 2008b] and 
phyllosilicates [Ehlmann et al., 2011] is spatially inhomoge¬ 
neous, with most fluvial features occurring at equatorial lat¬ 
itudes (north of 60° S). In addition, some equatorial regions 
such as Arabia and Noachis also exhibit low VN drainage 
densities, despite the fact that impact crater statistics indi¬ 
cate that the terrain is mainly Noachian [Tanaka, 1986]. 

In comparisons between climate models and the geologic 
evidence for early Mars, the focus is generally placed on the 
surface temperatures obtained, based on the argument that 
mean or transient temperatures above 0 °C are required to 
explain the observed erosion. Equally important, but less 
explored, is the extent to which the spatial distribution of 
fluvial features can be explained by climate models. The 
VN distribution has been comprehensively mapped in recent 
years [Hynek et al., 2010], so it constitutes a particularly 
important constraint on models of the early Martian water 
cycle. Here we begin by exploring a few of the most obvious 
explanations for why VNs are observed on some regions of 
Mars’ ancient crust but not on others. 

One seemingly obvious explanation for the predominance 
of VNs at equatorial latitudes is that early Mars was 
warmest at the equator, like Earth today. In an analo¬ 
gous climate to Earth’s with abundant water sources, this 
would increase either precipitation or the frequency of melt¬ 
ing events at the equator, and hence the erosion rate. How¬ 
ever, this explanation for the VN distribution is problematic, 
for two reasons. First, if early Mars had a thicker C0 2 at¬ 
mosphere of 1-2 bar [Phillips et al., 2001], equator-pole heat 
transport by the atmosphere would have been much more 
effective than it is today, which would have decreased the 
latitudinal temperature gradient. Second, Mars’ obliquity 
evolves chaotically due to secular orbital perturbations. At 
obliquity (f) = 41.8° (the most probable value predicted over 
timescales greater than a few Gy [Laskar et al., 2004]), the 
annual mean insolation for Mars 3.8 Gy at the poles is only 
30.0 W m -2 less than that at the equator, while the annual 
maximum (diurnal mean) insolation is 153.6 W m~ 2 greater 
(Figure 1). Both these effects argue against a warmer equa¬ 
torial climate as the reason for the observed valley network 
distribution. 

What about alternative explanations? Another possibil¬ 
ity is that differences in preservation dominate the spatial 
distribution, even on Noachian terrain. Difficulties in the de¬ 
tection of open-basin lakes poleward of 60° S has been noted 
[Fassett and Head, 2008b] and hypothesised to be a result of 
dust mantling at high southern latitudes [Soderblom et al., 
1973]. The presence of pedestal craters poleward of 40° S 
also suggests that modification of high latitude terrain by ice 
in the Amazonian may have been significant [Kadish et al., 
2009]. However, because the Martian obliquity varies chaoti¬ 
cally on geological timescales, even in the Amazonian water 
ice should have been transported from the poles to equa¬ 
torial regions and back again many times over [Madeleine 
et al., 2009]. If ice mantling has destroyed VNs that formed 
south of 60° S, therefore, it is not clear why it has not also 
destroyed them at equatorial latitudes. 

The third and most interesting possibility is that at least 
some of the VN spatial heterogeneity is a result of H 2 0 
supply limitations, due to regional differences in Noachian 
precipitation and/or sublimation/evaporation rates. If this 
is the case, it may provide clues as to the nature of the early 
Martian hydrological cycle, and hence the climate. To in¬ 
vestigate this idea, we have performed a range of 3D GCM 
simulations under both cold and warm conditions. 


3. Method 

To compare various climate scenarios without constraints 
based on the specific combination of greenhouse gases avail¬ 
able or the faintness of the young Sun, we assumed a C0 2 - 
dominated atmosphere as in Wordsworth et al. [2013] but 
studied the effects of a) progressively increasing the solar 
flux F from the value of 75% present-day predicted for 
~3.8 Ga from standard solar models and b) adding a gray in¬ 
frared absorption coefficient k. We emphasize that such sim¬ 
ulations are nominally unphysical: a significantly brighter 
young Sun is possible, but conflicts with both theoretical 
models of stellar evolution on the main sequence and obser¬ 
vations of nearby Sun-like stars (e.g., Minton and Malhotra 
[2007]). Intense long-term greenhouse warming by gaseous 
absorption or infrared cloud scattering in a dense C0 2 
atmosphere has been previously proposed [Pollack et al., 
1987; Forget and Pierrehumbert, 1997], but neither mech¬ 
anism currently appears physically viable [Kasting, 1991; 
Wordsworth et al., 2013; Forget et al., 2013]. Most recently, 
a H,,-C0 2 atmosphere to create a long-term warm and wet 
early Mars has been suggested [Ramirez et al., 2014], fol¬ 
lowing a similar proposal for H 2 -N 2 warming on the early 
Earth [Wordsworth and Pierrehumbert, 2013a]. This re¬ 
quires low H 2 escape rates combined with very high com¬ 
bined H., and CO a outgassing in the late Noachian, which 
appears petrologically difficult to achieve [Hirschmann and 
Withers, 2008]. Nonetheless, because the late Noachian geo¬ 
morphology is often taken as evidence for a warmer, wetter 
climate, performing ‘empirical’ 3D simulations of a warm, 
wet early Mars should still lead us to new insights. 

For simplicity, we focus on two end-member cases here. 
In the first “cold and icy” case, a solar flux of F — 0.75Fo = 
441.1 W m -2 is used as in Wordsworth et al. [2013] to match 
Mars’ received flux approximately 3.8 Ga [Gough, 1981] and 
no additional gray-gas atmospheric absorption is included. 
In the second “warm and wet” case, F is chosen to yield 
global mean temperatures sufficient for liquid H,0 precipita¬ 
tion (rain) in the southern highlands, and a northern ocean 
is assumed to be present. In both cases, the standard obliq¬ 
uity was taken to be 41.8° but values ranging from 10° — 55° 
were studied. Key model parameters are given in Tables 1 
and 2, and further details of the model setup are described 
below. 

3.1. Equilibrium temperature curves and valley 
network drainage density plot 

To produce the equilibrium temperature curves in Fig. 1, 
we assumed a constant planetary albedo A = 0.2 and calcu¬ 
lated net insolation F(L S , A) as a function of season angle 
L s and latitude A using standard expressions (e.g., Pierre¬ 
humbert [2011]) and neglecting the diurnal cycle. The valley 
network density plot was calculated using the Line Density 
tool in the geographical information system software pack¬ 
age ArcMap. Network density was calculated by considering 
a radius of 200 km around each pixel, summing the lengths 
of all the valley networks contained within it, and divid¬ 
ing this number by the area of the circle. Valley network 
outlines were taken from Hynek et al. [2010], and geologic 
units were taken from Scott and Tanaka [1986], Greeley and 
Guest [1987] and Tanaka and Scott [1987]. 

3.2. General model features 

To produce the 3D simulation results we used the LMD 
generic climate model [Wordsworth et al., 2010b]. Briefly, 
the model consists of a finite-difference enstrophy-conserving 
dynamical core that solves the meteorological primitive 
equations [Sadourny, 1975] coupled to physical parameter- 
izations for the planetary boundary layer and atmospheric 
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radiative transfer. The latter is solved using a two-stream 
approach [Toon et al., 1989] with optical data derived from 
line-by-line calculations [Rothman et al., 2009] and inter¬ 
polated collision-induced absorption data [Clough et al., 
1989; Gruszka and Borysow, 1998; Baranov et al., 2004; 
Wordsworth et al., 2010a]. We use the same spectral res¬ 
olution as in Wordsworth et al. [2013] (32 infrared bands 
and 36 visible bands) but an increased spatial resolution of 
64 x 48 x 18 (longitude x latitude x altitude). 

As before, the model contains three tracer species: C0 2 
ice, H,0 ice/liquid and H 2 O vapour. Local mean CO 2 and 
HoO cloud particle sizes were determined from the amount 
of condensed material and the number density of cloud con¬ 
densation nuclei [CCA], which was set at 10 5 kg' 1 for C0 2 
ice, 5 x 10 5 kg -1 for H 2 0 ice and 10 7 kg -1 for liquid H 2 0 
droplets. The value of [CCN] is extremely difficult to con¬ 
strain for paleoclimate applications. The values we have 
chosen for H a O are simply best estimates based on compar¬ 
ison with observed values on Earth [Hudson and Yum, 2002]. 
The dependence of climate on the assumed value of [CCN] 
for both C0 2 and H a O are studied in detail in Wordsworth 
et al. [2013] and Forget et al. [2013]. 

Convective processes and cloud microphysics were mod¬ 
eled in the same way as previously, with the excep¬ 
tion of H,0 precipitation. This was modeled using the 
parametrization of Boucher et al. [1995], with the termi¬ 
nal velocity of raindrops calculated assuming a dynamic at¬ 
mospheric viscosity and gravity appropriate for C0 2 and 
Mars, respectively. The Boucher et. al. scheme is more 
physically justified than the threshold approach we used 
in Wordsworth et al. [2013]; nonetheless in sensitivity tests 
(Figure 4b) we found that the two schemes produced com¬ 
parable time-averaged results. 

3.3. Warm and wet simulations 

In the warm and wet simulations, a 1 bar average sur¬ 
face pressure was assumed. Estimates of the maximum 
atmospheric CO a pressure during the Noachian are of or¬ 
der 1-2 bar [Phillips et al., 2001; Kite et al., 2014], with 
one study [Grott et al., 2011] pointing to lower values 
(~ 0.25 bar) based on outgassing models that incorporate 
the lower estimated oxygen fugacity of typical Martian mag¬ 
mas [Hirschmann and Withers, 2008]. 

The topography was adjusted so that the minimum al¬ 
titude was -2.54 km, corresponding to the putative north¬ 
ern ocean shoreline based on delta deposit locations from di 
Achille and Hynek [2010]. All regions at this altitude were 
then defined as ‘ocean’: surface albedo was set to 0.07 and 
an infinite water source was assumed (Figure 2). In addition, 
the thermal inertia was adjusted to 14500 tiu, correspond¬ 
ing to an ocean mixed layer depth of approximately 50 m 
[Pierrehumbert, 2011]. Horizontal ocean heat transport was 
neglected. Including it would have required a dynamic ocean 
simulator, which was outside the scope of this study. How¬ 
ever, at 1 bar atmospheric pressure the meridional atmo¬ 
spheric heat transport is already effective, suggesting that 
inclusion of ocean heat transport would not significantly al¬ 
ter the main conclusions. 

On land, runoff of liquid water was assumed to occur 
once the column density in a gridpoint exceeded q m ax,surf = 
150 kg m -2 , while the soil dryness threshold was taken to 
be 75 kg m -2 , corresponding to subsurface liquid water lay¬ 
ers of 15 cm and 7.5 cm, respectively. These values were 
chosen based on Earth GCM modeling [Manabe, 1969], and 
are reasonable first estimates given our simplified hydrol¬ 
ogy scheme and lack of knowledge of late Noachian subsur¬ 
face conditions. The sensitivity of the results to the runoff 
threshold are investigated in Section 4.4. 

The ice-albedo feedback was neglected in the warm and 
wet simulations by setting the albedo of snow/ice to the 
same value as that of rock (0.2). Its inclusion would have 


increased the chance of significant glaciation occurring in 
the simulations for a given solar flux or atmospheric infrared 
opacity and hence made a warm, wet Mars even more diffi¬ 
cult to achieve. An investigation of the stability of an open 
northern ocean on early Mars to runaway glaciation will be 
given in future work. 

Finally, for the idealised equatorial mountain simulations 
that were run to study the effect of Tharsis on the hydro- 
logical cycle, we used 

A) = 2 x 10V“ (x/40 ° )2 (1) 

with x 2 — ($ + IOO 0 ) 2 + A 2 , and longitude 9 and latitude 
A expressed in degrees. This describes a flat surface with a 
Gaussian height perturbation of full width at half maximum 
66.6° centered at 0° N and 100° W. 

3.4. Cold and icy simulations 

I 11 the cold and icy simulations, 0.6 bar average surface 
pressure was assumed. This value was chosen as the mini¬ 
mum necessary to ensure atmospheric stability against col¬ 
lapse of C0 2 on the surface for a wide range of obliquities. 
In some simulations, C0 2 ice appeared in specific regions on 
the surface all year round. However, secular increases in the 
total surface C0 2 ice volume after multiple years were not 
observed. 

The local stability of surface ice on a cold planet is con¬ 
trolled by two parameters: the accumulation rate and the 
sublimation rate. The sublimation rate can be calculated 
based only on the near-surface temperature, humidity, and 
wind speed. Calculating the accumulation rate, in contrast, 
requires a fully integrated (and computationally expensive) 
representation of the hydrological system. In Wordsworth 
et al. [2013], we handled this problem using an ice evolution 
algorithm. We ran the climate model with a full water cycle 
for several years, calculating the tendencies for sublimation 
and accumulation. We then multiplied them by a factor be¬ 
tween 10 and 100 and restarted the process until the system 
approached an equilibrium state. 

Here, the increased spatial resolution of our model makes 
even this method prohibitively expensive computationally. 
Instead, we calculated ice accumulation rates with a full 
water cycle in two scenarios with different ice initial condi¬ 
tions [ a) polar source regions and b) an H,0 source below 
-2.54 km (the ‘frozen ocean’ scenario)]. We also computed 
annual average potential sublimation rates for a wide range 
of parameters in simulations without a water cycle. 

We discovered that surface thermodynamics are far more 
important than the details of the large-scale circulation for 
governing the surface ice evolution in the model, leading to 
a situation where the spatial distribution of potential sub¬ 
limation rates is inversely correlated with that of snow ac¬ 
cumulation rates. This means that we can identify places 
where ice will stabilize in the long term far more rapidly 
than is possible using an ice evolution algorithm. Maps of 
potential sublimation can thus be used as a diagnostic for 
investigating the spatial distribution of stable ice. 

Quantitatively, we define the potential sublimation as the 
maximum possible H 2 0 mass flux from the surface to the 
atmosphere 

Spot = ~ Psat(T s ), (2) 

iXH 2 0-/ a 

assuming a) that the surface is an infinite source and b) 
that the atmospheric relative humidity is zero. Here Rh 2 o 
is the specific gas constant for H 2 0. The quantities |v|, 
T a and T s (the surface wind speed, the temperature of the 
lowest atmospheric layer and the surface temperature, re¬ 
spectively) were derived directly from the 3D model output 
and p S at (the saturation vapour pressure) was derived from 
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the Clausius-Clayperon relation. The drag coefficient Cd is 
defined (in neutrally stable conditions) as Cd = ( in[z/z 0 \ ) > 
where K. = 0.4 is the von Karman constant, z o is the rough¬ 
ness height, and 2 is the height of the first atmospheric layer. 
Cd was set to a constant 2.75 x 10~ 2 based on compari¬ 
son with the GCM value, which itself was calculated using 
a fixed roughness length of 10~ 2 m. Finally, we integrate 
Spot over one Martian year to give the theoretical maximum 
quantity of surface snow/ice that can be sublimated as a 
function of longitude and latitude. Our approach to calcu¬ 
lating Spot is accurate as long as the atmospheric thermody¬ 
namics is not strongly affected by the presence of H 2 0. This 
is the case for a cold early Mars, which has a surface moist 
convection number M of much less than unity [Wordsworth 
and Pierrehumbert, 2013b]. 

To produce estimates of transient melting in the cold, icy 
scenario, we performed simulations where we started from 
the baseline cold climate state and allowed the system to 
evolve after altering various parameters. The melting event 
simulations that included dust (Section 4.6) used the same 
method as in Forget et al. [2013], with the total atmospheric 
optical depth at the reference wavelength (0.7 /rm) set to 5. 
The value of e = 0.125 in the high eccentricity simulations 
was chosen with reference to the orbital evolution calcu¬ 
lations of Laskar et al. [2004]. The increase of solar flux 
from the baseline value of 0.75Fo to 0.8Fo in some of the 
simulations is justified given the uncertainty in the abso¬ 
lute timing of valley network formation as estimated from 
crater statistics [Fassett and Head, 2008a]. For some sim¬ 
ulations, the surface albedo was reduced. The rock albedo 
A r was taken to be that of basalt (0.1; plausibly the Mar¬ 
tian surface was richer in mafic minerals in the Noachian era 
[Mischna et al., 2013]), while the ice albedo A, was taken to 
be 0.3 (an approximate lower limit for ice contaminated by 
dust and volcanic ash deposition; Conway et al. [1996]). The 
simulations including S0 2 radiative forcing used the same 
correlated-fc method as used for the C0 2 -H 2 0 atmospheres, 
except that 10 ppmv of S0 2 was assumed to be present at 
all temperature, pressure and H 2 0 mixing ratio values cal¬ 
culated. S0 2 absorption spectra were obtained from the HI- 
TRAN database [Rothman et al., 2009], and the differences 
in line broadening due to the presence of C0 2 as the dom¬ 
inant background gas were neglected. Finally, the surface 
albedo was increased to A; and a thermal inertia appro¬ 
priate to solid ice (2000 tiu) was assumed in regions where 
more than 33 kg m -2 (3.3 cm) of surface H 2 0 was present 
[Le Treut and Li, 1991]. This implies a quite rapid increase 
of surface thermal inertia with snow deposition, which makes 
our results on the rate of summertime snowmelt in the cold 
scenario (Section 4.6) conservative. We also neglected fur¬ 
ther surface radiative effects such as insulation by a surface 
snow layer or solid-state greenhouse warming by ice, both 
of which could potentially cause some additional melting. 

4. Results 

4.1. Overview 

In our cold simulations, surface temperature is con¬ 
strained by a fully self-consistent climate model, whereas in 
the warm, wet simulations, we can choose surface tempera¬ 
ture by varying the solar flux and/or added gray opacity of 
the atmosphere. Fig. 3 shows annual mean surface temper¬ 
atures in the two standard scenarios, given obliquity 41.8°. 
In the cold case, surface pressure is 0.6 bar and solar flux 
is 0.75 times the present day value, while in the warm case 
the pressure is 1 bar and the solar flux is 1.3 times present 
day. Warm simulations where we used the best-estimate 
Noachian solar flux but added gray longwave opacity to the 
atmosphere showed broadly similar results (see Figure 9c). 

As can be seen, in both cases adiabatic cooling influences 
surface temperature. In the cold, icy scenario, the tem¬ 
peratures in the southern hemisphere are partly affected by 


seasonal condensation of C0 2 . In the warm, wet scenario, 
horizontal temperature variations in the northern and Hel¬ 
las basin oceans are low despite the absence of dynamical 
ocean heat transport, because heat transport by the 1 bar 
atmosphere is efficient. Global mean surface temperature 
is 225.5 K for the cold scenario and 282.9 K for the warm 
scenario. 

Fig. 4 shows VN drainage densities computed from high- 
resolution topographic, visible and infrared data [Hynek 
et al., 2010] alongside precipitation rates in the warm, wet 
scenario. For comparison, Fig. 5 shows snow/ice accumu¬ 
lation rates in the cold scenario. As can be seen, rainfall 
rates vary significantly with longitude as well as latitude. 
Over Noachian terrain, precipitation peaks occur in Arabia 
and around Hellas basin, including in the west and south¬ 
east where few VNs are observed. However, several regions 
with high VN drainage densities, such as Margaritifer Sinus, 
exhibit relatively low precipitation (approx. 10-20 times 
less than the peak value). In contrast, in the cold scenario 
(Fig. 5a-b), surface ice/snow mainly migrates to equatorial 
latitudes, where the majority of the VNs are observed. The 
snow deposition is inhomogeneous spatially, with peak val¬ 
ues of ~ 100 kg m~ 2 (Mars yr) _1 in the regions of highest 
elevation declining to under ~ 0.1 kg nr” 2 (Mars yr) -1 be¬ 
low 50°S. Snow accumulation rates in the cold, icy case are 
orders of magnitude lower than precipitation rates in the 
warm case due to the exponential dependence of the evapo¬ 
ration/sublimation rate on temperature. This of course also 
means that the equilibration time of the surface hydrological 
cycle is much longer when Mars is cold [Wordsworth et al., 
2013]. 

4.2. Warm and wet scenario: Latitude dependence 
of precipitation 

To understand the water cycle in our warm, wet simu¬ 
lations intuitively, we performed several idealized simula¬ 
tions. First, we calculated precipitation patterns in simula¬ 
tions with flat surface topography, infinite H,0 sources at 
every gridpoint, and entirely gray radiative transfer. Fig. 6 
shows the results of these ‘swamp planet’ simulations for a 
range of orbital obliquities. 

As can be seen, at all obliquities precipitation is highest 
near the equator in a narrow band of latitudes that varies 
with season. This band is the inter-tropical convergence 
zone (ITCZ). It is caused by regions of rapidly ascending air 
in the upward branches of the north- and southward Hadley 
cells [Holton and Hakim, 2013]. The second major feature 
of Fig. 6 is that for (j> = 25° or greater, significant precip¬ 
itation occurs at all latitudes, although in a pattern that 
appears distinct from the ITCZ. 

If instantaneous adjustment of the global mean flow to 
the solar forcing is assumed, moderate seasonal excursions 
of the ITCZ are possible [Lindzen and Hou, 1988]. This ex¬ 
plains the sinusoidal behavior of the precipitation at (j> = 10° 
and 25° in Fig. 6. However, when the latitude of peak solar 
forcing is too high, the double constraints of energy and an¬ 
gular momentum conservation cannot be met in the axisym- 
metric limit, and flow instability followed by eddy growth is 
expected in the winter hemisphere. In addition, the Hadley 
circulation in the summer hemisphere of seasonal axisym- 
metric flows becomes weak [Lindzen and Hou, 1988]. Given 
this, the ‘messier’ regions of high latitude precipitation in 
Fig. 6 at moderate to high obliquity must therefore be due 
to horizontally localized moist convection and/or eddy pro¬ 
cesses. Interestingly, precipitation maxima at high latitudes 
have also been modeled for the methane cycle on present-day 
Titan [Mitchell et al., 2006], although the physical processes 
at play in that case are somewhat different. 
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The localized nature of the high-latitude precipitation 
means that it is particularly sensitive to the presence of 
nearby water sources. Indeed, in simulations where we as¬ 
sumed Hellas Basin was dry, much less rainfall south of 60° S 
was observed (results not shown). Such simulations are 
likely physically inconsistent, however, because a large H 2 0 
reservoir combined with continual precipitation and ground- 
water infiltration should ensure the presence of liquid water 
in low geopotential regions across the planet’s surface [Clif¬ 
ford and Parker, 2001]. 


Note that in simulations we performed where Tharsis was 
removed, rainfall rates in the Margaritifer Sinus region were 
found to be high (results not shown). However, evidence 
from magnetic Held observations and geodynamics model¬ 
ing indicates that the bulk of Tharsis was already in place 
by the period of peak VN formation [Phillips et ah, 2001; 
Fassett and Head, 2011]. This suggests that erosion due to 
precipitation in a warm, wet climate is not the explanation 
for VN formation in Margaritifer Sinus. 


4.3. Warm and wet scenario: Tharsis rain shadow 

Another striking feature of the warm, wet simulations is 
the low precipitation around Margaritifer Sinus. Sensitivity 
studies indicated that this effect was caused by the presence 
of the Tharsis bulge. In our model, predominantly westerly 
surface winds are orographically lifted over Tharsis, causing 
enhanced rainfall on its western flank and reduced rainfall 
to the east. 

To understand the nature of this ‘rain shadow’ effect, 
we performed further calculations with simplified topogra¬ 
phy. The orientation of the effect was somewhat perplex¬ 
ing initially, because the annual mean equatorial surface 
wind direction was easterly in the flat simulations (as on 
Earth). To understand why the rain shadow appeared on the 
east side of Tharsis, we first performed 3D simulations us¬ 
ing the idealised Gaussian topography described previously. 
To decrease model complexity and gain greater insight, we 
assumed a dry, pure 1 bar C0 2 atmosphere, used non¬ 
scattering gray gas radiative transfer ( ki w = 2 x 10 4 , k 3W = 
0) and set obliquity to 0°. The choice of zero obliquity was 
made to remove the effects of seasonality from the results. 

Fig. 7 shows the results in terms of the vertical and hor¬ 
izontal velocity at a mid-tropospheric cr-level in the model 
(approximately 0.7 bar). As can be seen, the vertical veloc¬ 
ity distribution is zonally asymmetric at the equator, peak¬ 
ing on the west side of the topographic perturbation. This 
behaviour is qualitatively similar to that found in a simpler 
quasi-linear model of the present-day Martian atmosphere 
(see Webster [1977], Fig. 10). The zonal wind is westerly 
at most latitudes, which fits the picture of a predominately 
eastward rain shadow. 

Insight into a possible wave origin for the phenomenon 
can be gained by studying the classic Gill solution to the 
problem of equatorial waves forced by localised heating [Gill, 
1980]. Elevated topographic regions like Tharsis appear as 
heat sources to the atmosphere because of the same adia¬ 
batic cooling effect responsible for the icy highlands phe¬ 
nomenon in a cold climate: adiabatic cooling decreases the 
atmospheric temperature with altitude, but the solar energy 
input to the surface generally remains constant or increases 
with altitude. 

To get a sense of the linear response of the equatorial 
/3-plane system to simultaneous topographic and thermal 
forcing without a full eigenfunction analysis, in Fig. 8 we 
have plotted the u(x, ?/)-component of the (dimensionless) 
Gill solution for both Kelvin and Rossby modes 

u= ^[q 0 {x) + q2{x){y 2 -Z)}e~^ y2 (3) 


at the equator (y = 0), with go(x) and q 2 (x) defined as in 
equations (4.2) and (4.7) of Gill [1980]. We assume damping 
e = 0.05 and a width of forcing region L = 2. In Fig. 8b, 
the equatorial vertical velocity response w that would result 
if a topographic perturbation of the form H(x) = e~ k x 
were present [i.e., w ~ H'(x)u(x, 0) = —2k 2 xH(x)u(x,0) 
based on an assumption of incompressibility] is plotted. As 
can be seen, because of the zonal asymmetry of the Kelvin 
and Rossby wave responses, the vertical velocity can become 
strongly enhanced on the western side of the topographic 
perturbation. 


4.4. Warm and wet scenario: Sensitivity studies 

To test the robustness of our results, we also studied the 
sensitivity of our precipitation plots to variations in various 
poorly constrained model parameters. Results are shown in 
Figure 9. As can be seen, the same broad pattern is seen in 
all cases studied, with small differences on regional scales. 
Use of the simplified precipitation scheme from Wordsworth 
et al. [2013] (Fig. 9b) causes an increase in the global mean 
value, mainly due to increases in the north and over Tharsis. 
Keeping the solar flux at 0.75 Fo but adding gray longwave 
opacity to the atmosphere (Fig. 9c) decreases the global 
mean precipitation relative to the standard case, as expected 
from energetic considerations [Pierrehumbert, 2011]. It also 
causes precipitation to become more localized in specific re¬ 
gions. Changing the runoff threshold or number of cloud 
condensation nuclei (Fig. 9d-e) also has no major effect. 
Finally, decreasing obliquity to 25° slightly increases the 
amount of precipitation in the south and north and decreases 
it over Tharsis and Arabia Terra, but otherwise has little ef¬ 
fect (Fig. 9f). 

4.5. Cold and icy scenario: ice migration and annual 
potential sublimation 

Our high-resolution simulations yield broadly similar re¬ 
sults to those reported in Wordsworth et al. [2013]. Under 
higher atmospheric pressures than today, adiabatic cooling 
of the surface due to increased thermal coupling with the 
denser atmosphere causes surface temperatures to decrease 
with height on Mars [Forget et al., 2013; Wordsworth et al., 
2013]. This causes the equatorial highlands to be among the 
coldest regions of the Martian surface, leading to enhanced 
stability of snow deposits there (the ‘icy highlands’ effect 
[Wordsworth et al., 2013]). 

In Fig. 5, net surface H 2 0 snow/ice accumulation in 
simulations with a water cycle is plotted along annual po¬ 
tential sublimation in simulation with fixed atmospheric 
relative humidity but no active water cycle. As can be 
seen, at the most probable obliquity 3.8 Ga from Laskar 
et al. [2004] ( <j> = 41.8°), increased sublimation due to the 
high peak summertime temperatures causes the regions be¬ 
tween ~ 50° S and ~ 80° S to remain comparatively ice-free 
(Fig. 5a) and snow accumulation is spatially correlated with 
the VN drainage density (see also Fig. 13). Interestingly, 
when a frozen ocean is assumed present below -2.54 km 
(Fig. 5c), precipitation (snowfall) occurs over wider regions 
than in the standard cold case, including south of Hellas 
and Argyre Basins. Nonetheless, over multiple years the re¬ 
gion of peak snow/ice accumulation remains the southern 
equatorial highlands. 

As discussed in Section 3, the annual potential sublima¬ 
tion Spot (Fig. 5c) is closely correlated spatially with the 
snow accumulation rate. This demonstrates that surface 
thermodynamics is the key factor governing surface ice evo¬ 
lution in the model. The magnitude of S po t is significantly 
greater than the snow accumulation rate. This is expected 
because in reality the atmospheric relative humidity is non¬ 
zero and the surface is dry in many regions. S po t. indicates 
the regions where ice will stabilize in the long term, but it 
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gives only an approximate upper limit to the total rates of 
ice/snow transport. The close correlation of S po t with the 
snow accumulation rate in Figs. 5a-b nonetheless validates 
its use as a diagnostic for investigating the spatial distribu¬ 
tion of stable ice in other situations. 

In Figure 10, we show the variation of S po t with obliq¬ 
uity and surface pressure (Figure 10). In general, the effect 
of obliquity on surface ice stability lessens as the pressure 
increases, because the atmosphere becomes more effective 
at transporting heat across the surface (Figure lOf). At 
0.6 bar, the effect of obliquity is still fairly important, with 
lower values causing lower ice stability at the equator in gen¬ 
eral. At <j> = 25° (Figure 10b) S po t is lower (and hence the 
ice is more stable) in the equatorial highland regions than 
on most of the rest of the surface, although the main cold 
traps have become the south pole and peak of the Tharsis 
bulge. At <j> = 10° (Figure 10a), both the north and south 
poles are cold traps and water ice is driven away from equa¬ 
torial regions. Indeed, our simulations showed that at this 
obliquity C0 2 itself is unstable to collapse at the poles, even 
for a surface pressure as high as 0.6 bar. 

Transient melting of equatorial H,0 under a moderately 
dense or thin C0 2 atmosphere at low obliquity is there¬ 
fore only possible if ice migrates more slowly than obliquity 
varies. Variations in the Martian obliquity and eccentric¬ 
ity occur on 100,000 yr timescales [Laskar et ah, 2004]. Ice 
migration rates are difficult to constrain because of the non¬ 
linear dependence of sublimation rate on surface tempera¬ 
ture, but the mean values in our model for 0.6 bar C0 2 are 
of order 1 kg m -2 (Mars yr) -1 or ~ 1 mm (Mars yr) -1 , 
leading to ~10,000 Mars yr for the transport of 10 m global 
averaged equivalent of H 2 0. Unless the Noachian surface 
water inventory was orders of magnitude higher than the 
~ 34 m estimated to be present today, which may be un¬ 
likely [Carr and Head, 2014], the equator would therefore 
have been mainly dry at low obliquity. Hence in the cold, 
icy scenario for early Mars, equatorial melting events should 
only occur when the obliquity is 25° or greater and/or at¬ 
mospheric pressure is high. 

4.6. Cold and icy scenario: episodic melting events 

If early Mars was indeed mainly cold, some mechanisms 
must still have been responsible for episodic melting. Possi¬ 
ble candidates include seasonal and diurnal effects [Richard¬ 
son and Mischna, 2005; Wordsworth et al., 2013], positive ra¬ 
diative forcing from atmospheric dust and/or clouds [Forget 
and Pierrehumbert, 1997; Pierrehumbert and Erlick, 1998; 
Urata and Toon, 2013], orbital variations, meteorite impacts 
[Segura et al., 2008; Toon et ah, 2010] and S0 2 /H 2 S emis¬ 
sion from volcanism [Postawko and Kuhn, 1986; Mischna 
et al., 2013; Halevy and Head, 2014; Kerber et al., 2015]. 

Episodic impact events were an indisputable feature of 
the late Noachian environment, and they may also have had 
a significant effect on climate. Testing of the hypothesis that 
impacts caused melting sufficient to carve the valley net¬ 
works requires 3D modeling of climate across extreme tem¬ 
perature ranges, which we plan to address in future work. 
However, it is worth noting that the stable impact-induced 
runaway greenhouse atmospheres for early Mars proposed 
recently [Segura et al., 2012] are highly unlikely. The con¬ 
clusion of runaway bistability for early Mars is based on 
the assumption of radiative equilibrium in the low atmo¬ 
sphere, which requires unphysically high supersaturation of 
water vapour [Nakajima et al., 1992]. Under more realis¬ 
tic assumptions for atmospheric relative humidity, the small 
decrease in outgoing longwave radiation with surface tem¬ 
perature past the peak value in a clear sky runaway green¬ 
house atmosphere is not sufficient for hysteresis under early 
Martian conditions [Pierrehumbert, 2011]. 


We have systematically studied all the other effects men¬ 
tioned above (Fig. 11; see also Wordsworth et al. [2013]; 
Forget et al. [2013]; Kerber et al. [2015]). Alone, they are not 
sufficient to cause global warming sufficient to explain the 
observations. This is clear from Fig. lib, which shows the 
diurnal mean temperature vs. time averaged over the region 
0°W-120° W, 10°N-30° S (see box in Fig. 11a). The cyan 
line indicates the standard model, while the other colours 
indicate runs where effects were added separately [reduced 
surface albedo and increased atmospheric dust, blue; orbital 
eccentricity of 0.125 and solar constant of F = 0.8Fo, red; 
10 pprnv atmospheric S0 2 , green]. We find significantly 
lower warming by S0 2 than was reported in Halevy and 
Head [2014], primarily because our 3D model accounts for 
the horizontal transport of heat by the atmosphere away 
from equatorial regions (see also Kerber et al. [2015]). 

The maximum diurnal mean temperature for any of the 
studied cases is —20° C (for the increased eccentricity and 
F = 0.8-Fo run), implying very limited melting due to any 
single climate forcing mechanism. In combination, however, 
these effects can cause summertime mean temperatures to 
rise close to 0°C in the valley network regions. This is shown 
by the black line in Fig. lib; the gray envelope indicates 
the maximum and minimum diurnal temperatures. 

Figure 12 shows the average quantity of surface liquid wa¬ 
ter in the maximum forcing case as a function of time, as¬ 
suming a melting point of 0°C (black line) and -10°C (blue 
line). Uncertainty in the freezing point of early Martian ice 
due to the presence of solutes means that a melting threshold 
of-10°C is reasonable in a mainly cold scenario [Fairen et al., 
2009]. As can be seen peak melting reaches 5-7 kg m -2 in 
summertime in the briny case. This is still significantly be¬ 
low the 150 kg m -2 runoff threshold used in the warm and 
wet simulations, which itself was based on estimates from 
warm regions on Earth lacking subsurface permafrost. 

Hoke et al. [2011] estimate that intermittent runoff rates 
of a few cm/day on timescales of order 10 s to 10 7 years are 
necessary to explain the largest Noachian VNs. Our melting 
estimates are somewhat below these values, even if runoff 
is assumed to occur almost instantaneously after melting. 
Nonetheless, the exponential dependence of melting rate on 
temperature means that only a slight increase in forcing be¬ 
yond that used to produce Fig. 12 would bring the model 
predictions into the necessary regime. Furthermore, the de¬ 
tails of melting processes in marginally warm scenarios may 
still allow for significant fluvial erosion, as evidenced by ge- 
omorphic studies of the McMurdo Dry Valleys in Antarctica 
[Head and Marchant, 2014] and recent modeling of top-down 
melting from a late Noachian icesheet [Fastook and Head, 
2014]. Incorporation of a more sophisticated hydrological 
scheme in the model in future incorporating e.g. hyporheic 
processes will allow these issues to be investigated in more 
detail. 

In this analysis, we do not claim to have uniquely iden¬ 
tified the key warming processes in the late Noachian cli¬ 
mate. Nonetheless, we have demonstrated that melting in 
a cold, icy scenario is possible using combinations of phys¬ 
ically plausible mechanisms. In contrast, achieving contin¬ 
uous warm and wet conditions is much more difficult, as 
demonstrated by our simulations, which require an increased 
solar flux compared to present day (in conflict with basic 
stellar physics) and/or intense greenhouse warming in addi¬ 
tion to that provided by a 1-bar C0 2 atmosphere. 

4.7. Spatial correlation of valley networks with 
rainfall and snow accumulation rates 

Because the model predictions for precipitation in the 
warm, wet regime and snow accumulation in the cold, icy 
regime are quite different, it is interesting to compare the 
spatial correlations with the Hynek et al. valley network 
drainage density in each case. To do this, we first filtered 
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the data by excluding points on terrain that was not dated 
to the Noachian era. Polar areas and younger terrains were 
excluded from the analysis because absence of observable 
VNs does not necessarily indicate that they never formed in 
these regions. 

Figure 13 shows the resulting normalized logarithmic 
scatter plots of the VN drainage density (Fig. 4a), a) precip¬ 
itation in the standard warm and wet case (Fig. 4b) and b-c) 
snow accumulation in the cold and icy cases (Fig. 5a and 
b). For the warm and wet case, the correlation coefficient 
is -0.0033 and the p-value is 0.93, indicating the correla¬ 
tion is not statistically significant. For the two cold and icy 
cases, in contrast, the correlation coefficients are 0.23 and 
0.19 while the p -values are both < 0.001. 

The key reasons for the poor correlation in the warm, wet 
case appear to be a) the lack of precipitation around Margar- 
itifer Sinus and b) the high precipitation on terrain south of 
40°S, where few VNs are present. This is clear from Fig. 14, 
which shows the regions of (dis)agreement between the VN 
drainage density and precipitation/snow accumulation maps 
across the surface. Fig. 14 was produced by creating a mask 
from the VN data with a dividing contour defined such that 
50% of the surface was denoted as containing VNs. Similar 
masks were then created from the warm and wet precipita¬ 
tion and cold and icy snow accumulation data, and for each 
plot the masks were superimposed to highlight regions of 
agreement and disagreement. 

Figure 14 highlights the results discussed previously in 
Section 4.1: in our model the warm and wet scenario dis¬ 
agrees most with the VN distribution in Margaritifer Sinus 
and south of Hellas Basin. The cold and icy scenario, in 
contrast, has closer agreement around the equator but still 
some disagreement on terrain south of 40° S. In the cold and 
icy case, this discrepancy is somewhat artificial, because 
snow equilibration also varies significantly with obliquity. 
In the warm and wet case, the basic features of the pre¬ 
cipitation map are much less sensitive to obliquity changes 
(Sections 4.2-4.4). 

The lack of correlation between warm and wet precipi¬ 
tation and the VN distribution when compared to the cold 
and icy result is striking. Because prediction of precipitation 
patterns using GCMs is non-trivial even for the present-day 
Earth, further studies in future using a range of convection 
schemes and model dynamical cores will be important. If 
the dynamical effects we have analyzed here turn out to be 
independent of all model details, this clearly will constitute 
a significant piece of evidence against a warm and wet origin 
for many of the valley networks. 

5. Discussion 

We have presented the first direct comparison of warm, 
wet and cold, icy early Mars scenarios in a 3D climate model. 
Because of their generality, 3D models allow hypotheses for 
early Mars to be compared with the geological evidence and 
tested for internal consistency to a far greater extent than is 
possible with ID radiative-convective or 2D energy-balance 
models. 

In our warm and wet simulations, the precipitation pat¬ 
terns do not match the observed VN distribution closely. 
In particular, high precipitation is observed around Hellas 
and in Arabia Terra (where Noachian-era VNs are scarce). 
However, low precipitation is observed in Margaritifer Sinus 
(where the integrated VN drainage density is nonetheless 
high). 

In the cold and icy simulations, spatial variations in sup¬ 
ply rates are large, but the match between ice/snow accumu¬ 
lation and the VN distribution is closer. This is especially 
true at the most probable obliquity value predicted for the 


late Noachian, where statistically significant spatial correla¬ 
tion is found between the two datasets. At obliquities less 
than around 20° and pressures less than around 0.5 bar, 
both CO a and H 2 0 collapse at the poles, and the surface 
likely becomes too dry at equatorial latitudes for any signif¬ 
icant melting or runoff to occur. Formation of a large south 
polar H,0 cap at low obliquity but higher pressure was al¬ 
ready studied in Wordsworth et al. [2013], and may be linked 
to the late Noachian - early Hesperian-era Dorsa Argentea 
Formation [Head and Pratt, 2001; Kress and Head, 2015]. 

Our results have implications for NASA’s Curiosity mis¬ 
sion. Specifically, we find relatively high snow accumulation 
rates in the Aeolis quadrangle in the cold, icy case, demon¬ 
strating that an ice-free northern ocean is not required to 
supply H 2 0 to Gale crater. A mainly cold scenario for fluvial 
alteration in Gale crater is consistent with in-situ geochem¬ 
ical analyses of sedimentary material in the Yellowknife Bay 
formation [McLennan et al., 2014; Grotzinger et al., 2014], 

Achieving the necessary transient warming to cause the 
required erosion through physically plausible mechanisms is 
an ongoing challenge for the cold scenario, although we have 
demonstrated that some melting is possible due to combi¬ 
nations of volcanism, dust radiative forcing and orbital vari¬ 
ations. Impacts represent another major potential melting 
mechanism that by nature are also transient (e.g., Toon et al. 
[2010]). A warm, wet early Mars requires an anomalously 
luminous young Sun or intense greenhouse warming from 
an unknown source, and hence is less plausible from a pure 
climate physics perspective. 

Clearly, there are many ways in which this work could 
be extended. The parametrizations we use for cloud mi¬ 
crophysics and precipitation were chosen for their simplicity 
and independence from Earth-tuned parameters. It would 
be useful in future to compare them with more sophisticated 
schemes such as are currently used for Earth climate stud¬ 
ies (e.g., Khairoutdinov and Randall [2001]). In addition, it 
has been suggested recently that the two-stream approach 
to radiative transfer calculations (which we use) may be in¬ 
accurate for C0 2 cloud scattering [Kitzmann et al., 2013]. If 
this is the case, it would mean that our peak temperatures 
in Fig. 11 are slightly overestimated. It would of course also 
mean that warm, wet conditions in a C0 2 -dominated atmo¬ 
sphere are even harder to achieve. This issue is not critical 
to our results, but it should be settled in the future using a 
3D multiple-stream radiative transfer code. 

For the cold scenario, future modeling work must focus 
on the construction of self-consistent scenarios for melting 
episodes to explain the required erosion quantitatively. More 
sophisticated modeling of hydrology on a cold early Mars 
(e.g., Fastook and Head [2014]) will also lead to a better un¬ 
derstanding of the long-term evolution of water on the sur¬ 
face and beneath it. Further work should also be performed 
to reconcile the late Noachian-era equatorial geomorphology 
with the southern Hesperian-era Dorsa Argentea Formation 
in future. 

Finally, because our simulations make specific predictions 
for the spatial distribution of VNs and other fluvial features 
under different early climate scenarios, they strongly mo¬ 
tivate further geological investigations. In the future, we 
suggest that special attention should be paid to morpho¬ 
logical and geochemical analyses of Noachian terrain where 
predictions differ the most between warm, wet and cold, 
icy scenarios (e.g., in Arabia Terra and the Margaritifer Si¬ 
nus quadrangle). More detailed orbital or in-situ analysis 
of these regions will lead to further insight into early Mars’ 
global climate, and hence ultimately the question of whether 
there were ever long-term conditions that could have allowed 
a surface biosphere to flourish. 
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Table 1. Parameters used in the warm and wet climate simulations. Bold indicates the values for the ‘standard’ simulations. 


Surface pressure 

1 bar 

Solar flux 

441.1, 764.5 W m” 2 

Added atmospheric gray mass absorption coefficient 

0.0, 2.0x 10” 4 m 2 kg” 1 

Orbital obliquity 

25.0°, 41.8° 

Orbital eccentricity 

0.0 

Runoff threshold 

1, 150 kg m” 2 


Table 2. Parameters used in the cold and icy climate simulations. Bold indicates the values for the ‘standard’ simulations. 


Surface pressure 

0.125, 0.6, 2 bar 

Solar flux 

441.1 W m” 2 , 470.5 W m” 2 , 

Orbital obliquity 

10.0°, 25.0°, 41.8°, 55.0° 

Orbital eccentricity 

0.0, 0.125 

Atmospheric S0 2 

0.0, 10.0 ppmv 

Dust optical depth at 0.7 (im 

0.0, 5.0 




Figure 1. (left) Insolation vs. latitude for Earth (black) 
and early Mars with obliquity 25.0° (red) and 41.8° 
(blue). Solid and dashed lines denote annual mean and 
annual maximum (diurnal mean) temperatures, respec¬ 
tively. For Earth the obliquity is 23.4°. In both cases an 
albedo of 0.2 and orbital eccentricity of 0.0 is assumed, 
while for early Mars the solar constant is taken to be 0.75 
times the present-day value, (right) The corresponding 
annual mean (solid) and maximum (dashed) surface tem¬ 
peratures derived from dry early Mars 3D climate simu¬ 
lations with 1 bar surface pressure. 
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Figure 2. Plot of the surface water distribution in the 
warm, wet simulations. Brown and blue shading repre¬ 
sents land and ocean, respectively. 
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Figure 3. Plot of annual mean surface temperatures in 
a) the standard cold scenario and b) the standard warm, 
wet scenario. Obliquity is 41.8° in both cases. Global 
mean surface temperature is 225.5 K for a) and 282.9 K 
for b). 
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b) 

Figure 4. a) Surface VN drainage density (data de¬ 
rived from Hynek et al. [2010]), with shading of the back¬ 
ground corresponding to terrain age (light, Amazonian; 
medium, Hesperian; dark, Noachian). b) Annual precip¬ 
itation over 1 Martian year in the warm, wet simulation 
[global mean R avg = 0.4 m/(Mars yr)[. The solid black 
line indicates the ocean shoreline. 
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Figure 5. a) Surface snow/ice accumulation after 
5 years in the default cold 3D climate simulation (obliq¬ 
uity 41.8°; pressure 0.6 bar) given a) H a O sources at the 
poles; global mean A avg = 5.3 kg/m 2 and b) an H 2 0 
source and flat topography below -2.54 km; global mean 
value A aV g = 3.2 kg/m 2 (cold, wet scenario), c) Annual 
potential sublimation for a) and b). 
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Figure 6. Annual precipitation (snow and rain) in 
idealised early Mars simulations with flat topography, 
non-scattering gray gas radiative transfer (m w = 2 x 
10 -4 m 2 kg -1 , k 3W = 0.0 m 2 kg -1 ), uniform surface 
albedo of 0.2, solar flux F = 1.2Fo, mean atmospheric 
pressure of 1 bar, no C0 2 condensation, and an infinite 
surface H.,0 reservoir at every gridpoint. a) obliquity 
</>= 10°; b) 0 = 25°; c) <^ = 40°; d) <j> = 55°. 
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Figure 7. Annual mean vertical velocity (filled con¬ 
tours) and horizontal velocity (black arrows) at the 8th 
model level (approx, p = 0.7 p S urf) in an idealised dry, 
gray simulation with obliquity of 0° and topography rep¬ 
resented by equation described in the main text (black 
solid lines). 
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Figure 8. a) Zonal velocity in the Gill solution to the 
linear shallow water equations on an equatorial /3-plane 
with local, spatially symmetric heating . b) Vertical ve¬ 
locity implied from a) due to mass conservation given the 
presence of a Gaussian topographic perturbation. 
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Figure 9. Annual precipitation in warm, wet runs with 
ocean assumed as in Fig. 2 and standard parameters as in 
Table 1. a) Standard parameters as defined in the main 
text, b) a simplified threshold precipitation scheme as in 
Wordsworth et al. (2013), c) solar flux 0.75Fo but added 
gray gas atmospheric mass absorption coefficient Ki w = 
2 x 10 -4 m 2 /kg, d) runoff threshold set to 1 kg m -2 , e) 
variable cloud particle sizes with [CCN] set to 1 x 10 5 for 
all particles, f) obliquity (j> = 25°. Global mean values 
( Ravg ) are 0.41, 1.0, 0.31, 0.32, 0.65, 0.47 m/(Mars yr) 
for a-f), respectively. 
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Figure 10. Annual potential sublimation in cold simula¬ 
tions with a) mean atmospheric pressure Pco 2 = 0.6 bar, 
obliquity <j> = 10° and eccentricity e = 0, b), Pco 2 = 
0.6 bar, <j> = 25° and e = 0 c), Pco 2 = 0.6 bar, <j> = 41.8° 
and e = 0, d), pco 0 = 0.6 bar, (j> = 55° and e = 0 e), 
Pco, = 0.125 bar, <j) = 41.8° and e = 0, f) pco 2 = 2 bar, 
4> = 41.8° and e = 0 and g) pco 2 = 0.6 bar, <f> = 41.8° 
and e = 0.125. 
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Figure 11. Episodic warming in the cold, icy scenario, 
a) Annual maximum diurnal mean surface temperature in 
the standard cold simulation, b) Diurnal mean tempera¬ 
ture vs. time averaged over area denoted by the rectan¬ 
gle in the top panel for different forcing combinations [no 
forcing, cyan; reduced surface albedo and increased atmo¬ 
spheric dust, blue; orbital eccentricity of 0.125 and solar 
constant of F = 0.8Fo, red; 10 ppmv atmospheric S0 2 , 
green]. In black, a simulation with all of the previous ef¬ 
fects is shown for four years simulation time. There, the 
extremal temperatures are indicated by the gray shad¬ 
ing. The dashed line indicates the melting point of water, 
while the dotted line at -10° C indicates the point where 
up to 78% of briny water would remain liquid according 
to Fairen et al. [2009]. 
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Figure 12. Episodic melting in the cold, icy scenario. 
Surface liquid water amount vs. time averaged over the 
same region as in Fig. 11 for the most extreme forcing 
case (black line in Fig. 11). Here the black line shows 
surface liquid water in the case where 0° C is the melting 
point (fresh water), while the blue line shows the corre¬ 
sponding melt amount when -10° C is the melting point 
(brines). 
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Figure 13. Normalized logarithmic scatter plot of the 
VN drainage density from Hynek et al. [2010] vs. a) 
rainfall rate in the warm, wet scenario, b) snow accumu¬ 
lation in the cold, icy scenario with polar water sources 
and c) snow accumulation in the cold, icy scenario with 
low altitude water sources. From a-c), the correlation 
coefficients are -0.0033, 0.23 and 0.19, while the p -values 
are 0.93, < 0.001 and < 0.001, respectively. In all cases 
the total number of points is 644. 
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Figure 14. Comparison of the VN drainage density 
from Hynek et al. [2010] with a) rainfall in the warm, 
wet scenario, b) snow accumulation in the cold, icy sce¬ 
nario with polar water sources and c) snow accumula¬ 
tion in the cold, icy scenario with low altitude water 
sources. Green indicates regions of agreement between 
VNs and rain/snow. Orange indicates significant pres¬ 
ence of VNs but low rain/snowfall. Magenta indicates 
significant rain/snow but few VNs. Gray indicates low 
VN density and low rain/snowfall. Regions in white are 
where there was no data or the terrain is not Noachian 
age. 









